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We present the results of a systematic study using the density functional theory (within the local 
density approximation) of the effects of biaxial strain and composition on the self-diffusion of Si 
and Ge in Sii-xGe^ alloys diffusing by a vacancy mechanism. The biaxial strain dependence of the 
vacancy formation energy was reconfirmed with previous calculations. The effect of biaxial strain 
on the interaction potential energy between a substitutional Ge atom and a vacancy was calculated. 
The interaction potential energy included not only the ground state energies of the vacancy at 
different coordination sites from the Ge atom but also the migration energy barriers to jump from 
one coordination site to the adjacent. These calculations were used to estimate the change in 
the activation energy (due to biaxial strain) for the self-diffusion of Si and Ge in Si by a vacancy 
mechanism. The composition dependence of the vacancy formation energy was calculated. A 
database of ab initio migration energy barriers for vacancy migration in different local environments 
was systematically developed by considering the effect of the first nearest neighbor sites explicitly 
and the effect of the other sites by a mean field approximation. A kinetic Monte Carlo simulation 
based on the migration energy barrier database was performed to determine the dependence (on 
the composition) of the activation energy for the diffusion of Si and Ge in Sii-^Gex. A detailed 
study of the variation of the correlation factor with composition and temperature in Sii-^Gea; was 
performed using the results of the KMC simulation. These analyses constitute essential building 
blocks to understand the mechanism of vacancy mediated diffusion processes at the microscopic 
level. 



I. INTRODUCTION 

Silicon germanium technology is becoming increasingly 
popular in high frequency, low power applications p the 
principal reasons being the advancement in precision 
growth technologies^ and the compatibility of Sii_ x Ge x 
with the Si manufacturing processes along with such 
properties of Sii^Ge^ as the composition dependence of 
the band gap, the strain dependence of the carrier mobil- 
ity, and the increased dopant solubility in Sii-^Ge^ com- 
pared to Si. The abrupt change in the Ge concentration 
between the Si and the Sii-^Ge^ layers being a functional 
necessity in these devices, a key materials issue is that of 
interdiffusion in these layers. There have bee 
studies of this stress-coupled interdiffusion] 
Yet, a lot about the actual microscopic mechanisms be- 
hind these phenomena remains to be understood. Sim- 
ilarly the microscopic mechanisms responsible for the 
growth, composition, and the shape of Sii-aGe^ islands 
on Si, which find applications in areas like Si-based quan- 
tum dots, are also not well understood. From a techno- 
logical standpoint, understanding diffusion in Sii_ x Ge x 
is therefore important. From a scientific standpoint, the 
silicon germanium system presents an ideal and clean sys- 
tem (without the complications introduced by charged 
defect states) to further the theoretical understanding of 
interdiffusion in-random alloys in general. Quoting from 
a recent paper ,EJ "Theoretical treatments of self-diffusion 
in SiGe are uncharted areas - and the effect of strain even 



more so." These reasons have motivated us to perform 
a systematic and detailed first principles study of the 
Sii-^Ge^ system. 

In their recent paper, Zangenberg et aZ0 have pre- 
sented their results of a systematic experimental study 
of the variation of Ge self-diffusion in mono crystalline 
Sii-zGea; epi-layers as a function of composition (x) and 
biaxial strain. p(A similar study has also been reported by 
Strohm et aLlI3) These works represent an advancement 
over that presented by McVay and DuCharmeEJ in 1974, 
which studied the composition dependent Ge diffusion in 
polycrystalline Sii-^Gea;. These results have been used 
in the past as an input to empirically explain other exper^. 
imentally observed phenomena. For example, Baribeaucl 
used the Ge dependent diffusivity from Ref. [l3| to nu- 
merically solve the one dimensional Fick's diffusion equa- 
tion and compared the results with those of the expe: 
imentally determined ones. Similarly, Aubertine et al 
have used the Ge dependent diffusivity from Ref. |ll| in 
a commercial numerical solver to perform a similar com- 
parison with their experiments to provide an empirical 
explanation to the experimentally observed time depen- 
dent interdiffusivity in Si/Sii-zCe^ multilayers. Thus, 
although these results (of Refs. [Il],|l^,[l3|) have been valu- 
able in providing empirical insights into other phenom- 
ena, the reasons for these behaviors themselves have not 
been queried into from a fundamental level, to the best 
of our knowledge. We ae€ aware, however, of a recent 
paper by Venezuela et al£-& that has made an attempt in 
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this direction. 

Our present work is a part of a project intended 
to develop a fundamental understanding at the micro- 
scopic level ultimately, of the coupled strain relaxation 
and interdiffusion phenomenon in Si/Sii^Ge^. multi- 
layers. Towards this end, a fundamental understand- 
ing of the strain and composition dependent diffusivity 
in Sii^G&q. a s o bserved by the experiments mentioned 
previouslytl'tHlEJ is an essential prerequisite. Diffusion 
in SiGe has been postulated!!! to be mediated by point 
defects: vacancies, interstitials and by a point defect free 
mechanism - the concerted exchange mechanism— Very 
recently, another point defect, which the authorst£l have 
termed the fourfold coordinated defect (FFCD), has been 
suggested which could also be responsible for diffusion in 
SiGe. From their experimental observations, Fahey et 
al.ll suggest that at 1050 °C, the vacancy mechanism 
probably contributes to 60% - 70% of the Ge diffusion in 
Si, the rest being due to the interstitial mechanism. We 
note that the results presented in Ref. |l7| are for the dif- 
fusion of Ge in pure Si and so the relative contributions 
of the different mechanisms could be different for systems 
with different Ge concentrations. However, because va- 
cancies are among the important contributors, they are 
the focus of this article. 

Density functional theory (DFT) calculations have 
played a significant role in computational physics dur- 
ing the past few decades since the theory's formal incep- 
tion in the mid 1960s. The unknown nature of the ex- 
change correlation functional and the inability to make 
progressively more accurate approximations to the same 
(as would be possible, for example, in a wave function 
based method), however, has been one of the main is- 
sues concerning the practical application of the DFT. 
The search for better exchange correlation functionals is 
an active area of research in the theoretical physics com- 
munity. The popular local density approximation (LDA) 
and the more computationally expensive (but not neces- 
sarily more accurate) generalized gradient approximation 
(GGA) , unfortunately, have been unable to reproduce ex- 
perimentally observed values of activation energy of diffu- 
sion in Si, the discrepancytfl Ueiug as high as leV. Quan- 
tum Monte Carlo techniques ,E2lEj which circumvent the 
problem due to the exchange correlation functional, are 
gaining popularity . However, because the LDA based 
DFT is definitely one of the most advanced computa- 
tional tools available for the systems of the size that we 
would like to study, we have used it in this present study. 
Because we are aware of this discrepancy between the 
theoretical prediction and the experimental values, and 
because we have restricted this present article to only 
the vacancy mechanism, we prefer to refrain from mak- 
ing very strong comparisons of our results with those of 
experiments, leaving such comparisons to the future un- 
til we have resolved these outstanding issues. In spite 
of these limitations, we believe that our contribution is 
significant for the following two reasons: (i) Our results 
can be viewed as that of the behavior of a hypotheti- 



cal random binary alloy system (with the energetics pro- 
vided by the LDA) diffusing by a vacancy mechanism. 
This is of basic scientific interest, (ii) The infrastruc- 
ture that we have developed in this present work can 
be reused with little effort once accurate energetics be- 
comes available. This, along with similar analyses for 
other diffusion mechanisms can then be used to directly 
compare/predict experimental observations. 

Because diffusion is a thermally activated process, dif- 
fusivity can be characterized by a temperature indepen- 
dent term (the pre-exponential factor, Dq) and a tem- 
perature dependent term (the exponential of the activa- 
tion energy, E a , i.e., exp[— Ea/ksT] where ks and T are 
respectively the Boltzmann's constant and the temper- 
ature). For a defect mediated diffusion mechanism, the 
activation energy (E a ) is composed of the defect forma- 
tion energy and the rest, which we have termed as the 
activation-minus-formation (AMF) energy. We note that 
this has traditionally been called as the migration energy. 
The reason we have not called it the migration energy 
needs some explanation. We first consider the case of a 
tracer self diffusion. Shown in Fig. [j], for the sake of il- 
lustration, is the motion of a tracer in a two dimensional 
hexagonal lattice from an initial state to the final through 
the saddle point. The migration energy in this case has a 
very direct physical significance, namely, the energy dif- 
ference between the saddle point and the initial (ground) 
state. But when one considers anything other than a 
unary system, for example the diffusion of a tracer Ge in 
Si, one is unable to make a physically appealing corre- 
spondence to a migration process as in the case of a tracer 
self diffusion in a unary system. Specifically, considering 
the illustration shown in Fig. ^, for the Ge atom (filled 
circle) to be effectively displaced from its current position 
(labelled 1), the vacancy (filled square) has to move from 
its current position (labelled 2) to atleast the third coor- 
dination site from Ge (labelled 8) and return by another 
path (for example through sites labelled 6 and 3). As 
shown in Fig. ^, there are different energy barriers to get 
to different configurations: the barrier for the vacancy to 
get from the first coordination site (from Ge) to the sec- 
ond is different from that to get from the second to the 
third. The contribution to the activation energy (other 
than the effective defect formation energy) is not only 
from the complex collective action of all these different 
migration barriers but also, in this particular pair diffu- 
sion model, due to the energy required for the vacancy 
to get to the third coordination site, for example, so that 
it can return to the Ge atom from a different direction 
thereby causing a net motion of the Ge atom. These 
complex collective actions manifest in different forms, for 
example, as a temperature dependent correlation factor. 
Therefore, we felt the need to make the distinction from 
the term: migration energy. The AMF energy equals 
the migration energy in the microscopic sense only for 
the case of the tracer self^diffusion in a unary system. 
Although previous reportsEJ'EBE^I have suggested differ- 
ent measures of an effective migration barrier for these 
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FIG. 1: Migration energy for the motion of a tracer (shown 
as the target symbol) in a unary system has a direct physical 
correspondence: difference between the saddle point and the 
ground state energies. 




FIG. 2: Schematic of the Si structure. Open circles - Si; filled 
circle - Ge; filled square - vacancy. 



cases, we remark that we do not find any of them physi- 
cally enlightening. 

In order to study the effect of biaxial strain and 
composition on the diffusivity, one needs to study 
their effect on these three parameters, viz., the pre- 
exponential factor (Do), the vacancy formation energy, 
and the AMF eneca 



.Several previous first princi- 



ples calculationsEjCJu3'E£l have reported on the strain 
dependence of vacancy formation energy. In this arti- 
cle, we present our results where we have reconfirmed 
the biaxial-strain dependence of the vacancy forma- 
tion energv~.-Though we have seen general theoretical 
treatmentsEZrEa of the effect of strain on defect diffusion, 
we have not come across any reference reporting on the 
ab initio based calculation of the variation of the AMF 
energy (or the effective migration energy as it is gener- 
ally known) as a function of strain. In this article, we 
have calculated the biaxial strain dependence of the in- 
teraction potential energy between a substitutional Ge 
atom and a vacancy. The interaction potential energy 
included not only the ground state energies of the va- 
cancy at different coordination sites from the Ge atom 
but also the migration energy barriers for jumps from 
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FIG. 3: Interaction potential energy (in eV) between a substi- 
tutional Ge atom and a vacancy as a function of vacancy po- 
sition in relaxed Si from LDA calculations. (Lines are drawn 
as a guide to the eye.) 



one coordination site to the adjacent. We have used these 
calculations to estimate the change in the activation en- 
ergy (due to biaxial strain) for the self-diffusion of Si and 
Ge in Si by a vacancy mechanism. We then present our 
calculations of the Ge concentration dependence of the 
vacancy formation energy where we have used the clas- 
sic Boltzmann factor enhancement of the probabilities. 
We note that Venezuela et a/.EJ have used an approach 
similar to ours in thcir—tccciit paper. Earlier theoretical 
or numerical studi e d 2 jEjEjE 3 have only reported on the 
analyses of the concentration dependence of diffusivity 
in the low concentration regime, typically those corre- 
sponding to dopant concentrations. We have not found 
theoretical or numerical treatments of the concentration 
effects on diffusivity for higher solute concentrations like 
those found in Sii_ x Ge x alloy systems. We present the 
variation of the AMF energy as a function of Ge con- 
centration, which we have obtained from kinetic Monte 
Carlo (KMC) simulations using a migration energy bar- 
rier database calculated from first principles. The KMC 
simulations also enabled us to study the variation of the 
correlation factor as a function of Ge concentration and 
temperature. Such a study provides useful insight into 
the vacancy mediated diffusion mechanism in a random 
binary alloy arranged in a tetrahedral geometry. 

This article is organized as follows: In section |n[ we 
present the theoretical details of our calculations and the 
computational details of our simulations. In section III, 
we discuss our main results. In section IV, we conclude 



the article with a summary of the present work and also 
comment on the limitations of this work. 
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II. THEORETICAL AND COMPUTATIONAL 
DETAILS 

In this section, we present the theory behind our cal- 
culations and the main computational details of our sim- 
ulations. 



A. First principles calculations 

Our first principles calculations were performed 
using . the , .-plane- wave ultrasoft pseudopotential code 
VASPl23c3'Lj'Ea within the local density approximation 
(LDA). A 64- atom supercell with a kinetic,energy cut-off 
of 10 Hartree, and a 2 3 Monkhorst-PackES k-point sam- 
pling was used. Electronic energy convergence of up to 
2.7 x 10~ 5 eV was used and the structures were relaxed 
until the maximum force on any atom was less than 0.015 
eV/A. Saddle point configurations were determined us- 
ing the nudged elastic band (NEB) methodHZl Optimized 
Si and Ge lattice constants (computed by fitting the-tp* 
tal energy vs. the supercell volume to Murnaghan'sr 8 lr 9 l 
equation of state) were found to be 5.39 A and 5.62 A 
respectively. The vacancy formation energy and the va- 
cancy formation volume in Si (Ge) were found to be 3.31 
(1.88) eV, and -0.059 (-0.195)0 respectively where fl rep- 
resents the volume of a Si (Ge) atom. (We recall that the 
vacancy formation volume is the sum of the relaxation 
volume (Si: -20.73 A 3 ; Ge: -26.56 A 3 ) and the atomic 
volume (Si: 19.57 A 3 ; Ge: 22.23 A 3 ).) These values 
are expectecfexoiimarable to other recent first principle 
calculations .OOc3c3 Because of the low Si- vacancy for- 
mation volume_aiid because of Sii-^Ge^ being a model 
random alloy,OC3 the lattice parameters for Sii_ x Ge x 
were chosen by a simple rule of mixtures. Lattice param- 
eters so chosen were assumed to correspond to a strain 
relaxed state. Unless otherwise mentioned, all our calcu- 
lations were done in such a strain relaxed state. 



B. Kinetic Monte Carlo simulations 

The diamond lattice was generated in the computer 
memory by a four-dimensional integer array. The first 
three indices were used to reference the location (i.e., 
X,Y,Z "coordinates" ) of the cubic unit cell. The last in- 
dex was used to reference the particular atom among the 
eight in the unit cell referred to by the first-three indices. 
Because Sii-^Ge^ forms a random alloy,r 4 lr 5 l each mem- 
ber of this array was randomly designated as either a Si 
atom or a Ge atom and one randomly chosen member 
was designated as the vacancy. The relative numbers of 
the Si and Ge atoms were so chosen that the required 
composition was obtained. The displacement and the 
number of jumps performed by each of these atoms were 
recorded through out the simulation. Periodic bound- 
ary conditions were used so that an atom hopping out 



of one side reenters the system from the opposite side, 
essentially simulating an infinite system. 

Each KMC move consisted of the following five steps: 
(i) Obtaining the rates ri for the possible final configura- 
tions starting from the current configuration as the initial 
configuration, (ii) Generating a pseudo-random number 
7 € (0, 1]. (iii) Advancing the clockcj by — hi(7)/ ^ r i- 
(iv) Reconfiguring the system into one of the final con- 
figurations based on the random number generated in 
step (ii). (v) Updating the displacement and the num- 
ber of jumps of the vacancy and the atoms that have 
moved. We refer the reader to the original KMC paper 
by Voteio for details on the theory of the kinetic Monte 
Carlo algorithm. 

The supercell comprised of 50 x 50 x 50 cubic unit 
cells each containing eight lattice sites making up one 
million lattice sites. Random alloys of Sii-zCe^ with 
the concentration of Ge (x) varying from to 1 were 
used to study the effect of composition. A single va- 
cancy was used. (We note that the presence of one va- 
cancy in a 50 x 50 x 50 super cell, which we are forced 
to use due to computational limitations, results in an ex- 
tremely high concentration of vacancy compared to the 
real Sii-zGe^ system.) Three different random distribu- 
tions of Ge atoms were used for each composition and 
three different random number sequences were used for 
each distribution, thus making up nine samples for each 
composition. The results were averaged over all the nine 
samples. A billion vacancy hops were performed for each 
case. The scatter in the results among these nine samples 
was found to be extremely low. 



C. Effective vacancy formation energy calculation 

In this subsection, we explain how we use the classic 
Boltzmann factor to calculate the effective vacancy for- 
mation energy in the Sii^Ge^ alloy (where x denotes the 
atomic fraction of Ge). The strength of the interaction 
between a vacancy and a Ge atom is expectedly depen- 
dent on their relative positions. The interaction of the 
vacancy with the Ge atoms which are first nearest neigh- 
bors to the vacancy is stronger than with those which are 
second nearest neighbors. This second nearest neighbor 
interaction in turn is stronger than between those which 
are further away. We have therefore chosen three differ- 
ent forms to represent these three different interactions. 
For the strongest interaction, we define a function F. We 
denote as F(b) the drop in energy of the system when a 
vacancy is surrounded by b Ge atoms at the first nearest 
neighbor position to the vacancy. (For the Si structure, 
b, of course, ranges from zero through four.) For the in- 
teraction between the vacancy and Ge atoms that are at 
the second nearest neighbor positions to the vacancy, we 
use a linear expression for the drop in the energy with 
the number of Ge atoms in the second nearest neigh- 
bor position. We denote the constant of proportionality 
i.e., the drop in energy of the system for each Ge atom 
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in the second nearest neighbor position as S. Because 
the interaction between the vacancy and Ge atoms that 
are present beyond the second nearest neighbor positions 
is comparatively weak, we consider their effect through 
a mean field correction factor: M. M is the difference 
between the energy of a system with a vacancy whose 
first and second nearest neighboring positions are occu- 
pied by Ge atoms and all other positions are occupied by 
Si atoms and that of a system with a vacancy with Ge 
atoms in all the positions. We denote as E(n, b, x) the 
vacancy formation energy in a Sii-^Ge^ system (with a 
Ge concentration of x). Here, n denotes the total num- 
ber of Ge atoms in the first and second nearest neighbor 
positions to the vacancy and b denotes the number of Ge 
atoms that are in the first nearest neighbor positions to 
the vacancy. We then obtain the following expression for 
E(n, b, x) in terms of the vacancy formation energy in 

J Vf> 

above: 



pure Si (Ey) and the quantities F, S, and M defined 



E(n, b, x) = Epj ~ F(b) ~ (n - b)S 



Mx 



(1) 



If the distribution of Ge atoms is unaffected by vacan- 
cies, then, the probability p(n, b, x) of a vacancy being 
surrounded by n Ge atoms, b of which are first nearest 
neighbors to the vacancy and the rest [n — b) are second 
nearest neighbors is calculated in the following manner 
using the binomial Bernoulli distribution from elemen- 
tary probability theory: (Note: There are 4 first nearest 
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neighbor sites and 12 second nearest neighbor sites in 
the diamond lattice.) The probability of a vacancy being 
surrounded by b Ge atoms in the first nearest neighbor 
position is {^jx b (l — x) 4 ~ b (where x is the concentration 
of Ge). The probability of a vacancy being surrounded by 
(n — b) Ge atoms in the second nearest neighbor position 
is Ql h )x n - b (l - a;) 12 -( n - b ). The required probability 
p(n, 6, x) is therefore the product of the above two which 
simplifies to the following expression: 



P(n, b, x) 



12 



x n (l 



(2) 



The interaction between the Ge and the vacancies, 
however, affects their distribution. The probability 
p(n,b,x,T), which takes this interaction into account, 
is obtained by multiplying p(n, b, x) by the Boltzmann 
factor corresponding to the energy drop because of this 
interaction. We thus obtain the following expression for 
p(n,b,x,T): 

p(n, b, x, T) = p(n, b, x) exp[(F(6) + {n- b)S)/k B T] (3) 

where ks is the Boltzmann's constant and T is the tem- 
perature. We then express the effective vacancy forma- 
tion energy, (Ef(x, T)), as an average of the vacancy for- 
mation energies in the different environments weighted 
by their corresponding probabilities: 



(E f (x,T)) = 



4 / n 



16/4 



[J2p{n,b,x,T)xE(n,b,x)) +J2 (J^pM^T) x E(n,b,i 

I — 



.n=0 \6=0 



n=5 \b=0 



(4) 



where z is like the partition function: 



D. Theoretical calculation of correlation factor 
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£ 5>(n, b,x,T))+J2[ 5>(n, b, x, T) (5) 



7i=0 \6=0 



\b=Q 



We make the following two clarifications: (i)We have not 
included the mean field correction term M in the expres- 
sion for the Boltzmann factor in Eq. (||) because, M being 
independent of n or b, is factored out of the numerator 
and the denominator (z) in the expression for (Ef(x,T)) 
(Eq.( Eh) even if it is included, (ii) We have two terms 
in the RHS of Eqs. (Q) and (||) for the following simple 
reason: The number of first nearest neighbor Ge atoms 
(b) can be at most equal to the total number of Ge atoms 
in the first and the second nearest neighbor positions (n) 
when n is less than or equal to four. This is the first term 
on the RHS. The variable b can be at most equal to four 
when n is greater than four. This is the second term. 



The correlation factor (/) is definedefl as the ratio of 
the actual diffusion coefficient to the uncorrelated diffu- 
sion coefficient under the assumption that all the jumps 
are statistically independent of one another. The correla- 
tion factor provides a lot of insight into the microscopic 
mechanism of diffusion. In this sub-section, we give a 
brief outline of how the correlation factor is computed. 
We refer the reader to Ref. EJ| for further details. The 
correlation factor for the diffusion of a single impurity 
atom in a cubic crystal by the vacancy mechanism can 
be calculated using the expression: 



/ = 



1 



1 - (cos( 



(6) 



Here, (cos 8), which denotes the average of the cosine of 
the angle between successive impurity jumps, can be eval- 
uated as Tj cos Oj . Tj is the probability for the impurity 
to jump to the j th configuration. Oj is the angle formed 
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between the current impurity jump direction and the im- 
purity jump direction leading to the j th configuration. 
There is the implicit sum over the repeated index j. In 
the case of the diamond structure, j ranges over the four 
first nearest neighbors; i.e., the j th configuration results 
when the impurity jumps to the j th first nearest neighbor. 
Referring to Fig. 0, where the destination configurations 
(i.e., first nearest neighbors to impurity) are denoted by 
the numbers 2 through 5, the probabilities T 2 through T 5 
have been calculated by including jump sequences up to 
five vacancy hopsCJ in the following manner: We denote 
respectively by Vi, Vhi Vf, and vb, the jump rates for the 
following vacancy jump processes: (i) Vacancy exchanges 
positions with the impurity atom, (ii) Vacancy exchanges 
positions with the host atom without either breaking or 
forming a bond with the impurity atom, (iii) Vacancy ex- 
changes positions with the host atom and in the process 
forms a bond with the impurity atom, (iv) Vacancy ex- 
changes positions with the host atom and in the process 
breaks a bond with the impurity atom. (We explain the 
procedure for obtaining the values of the various jump 
rate s (vi, vh, vp, and vb) from our LDA calculations in 
Sec. [II A .) We denote as R(j,k) the probability for the 
impurity to jump to the j th first nearest neighbor posi- 
tion as a result of the vacancy performing k hops. (For 
example (referring to Fig. |^), one of the ways in which 
the impurity atom can jump to the first nearest neighbor 
position denoted as 2 as a result of the vacancy perform- 
ing three hops would be the following jump sequence of 
the vacancy: 2 to 7 followed by 7 to 2 followed by 2 to 1.) 
Using elementary probability theory, the various R(j, fe)'s 
can be computed to obtain: 



#(2,1) 
#(2,3) 
#(2,5) 



Vl + 3v B 



3 x 
9 x 



Vp 



vi + 3v B 
vi + 3v B 



v F + 3v H 



X 



v F -+ 

Vl 



vi 
1 

X - 

3v H 4 



3v B 



vp + 3v H vi + 3v B 
#(3,5) = #(4,5) = #(5,5) = 

vb vh 



(7) 
(8) 

(9) 



2 x 



vi + 3v B 

VF 



Vl 



X - 

3v H 4 



vp + 3vh vi + 3vb 



(10) 



Tj is then 



The R(j, fe)'s not listed above are all zero 
obtained by .simply summing R(j, k) over k from one 
through five.CJ We obtain the following expressions, in 
terms of the various R(j,k)'s, for the probabilities T 2 
through T5: 



T 2 
T 3 

Also, from Fig 
cos #5 = 1/3. 



#(2, 1) + #(2, 3) + #(2, 5) 
T 4 = T 5 = #(3, 5) 

j, cos6> 2 = —1 and cos6* 3 = cos 64 = 



(11) 
(12) 



E. Calculation of the correlation factor from the 
KMC simulation results 

The procedur e for calculating the correlation factor 
outlined in Sec. [I D is valid only for a single impurity 
atom migrating by a vacancy mechanism. Certain sym- 
metry requirements, which were implicitly used in the 
formulae presented there, are violated at higher impu- 
rity concentrations. In this subsection, we explain the 
procedure for calculating the correlation factor that is 
valid for any impurity composition, as long as there are 
a sufficient number of atoms to obtain a good statistical 
average. This procedure is a straightforward interpreta- 
tion of the definition of the correlation factor as applied 
to the results from the KMC simulation. 

From the definition of the correlation factor as the ratio 
of the actual diffusivity to the uncorrelated diffusivity 
and from the definition of the diffusivity as the ratio of 
the mean squared displacement (A 2 ) to 6r, where r is 
the time taken for the motion in the limit as r tends 
to zero, we obtain the correlation factor to be the ratio 
of the actual mean squared displacement to the mean 
squared displacement when the motion is uncorrelated. 
Symbolically, 



(X 2 



ctual 



(A 2 



random 



(13) 



From the random walk analysis, (A 2 ) ran dom = (A)A 2 , 
where (N) is the mean number of jumps and A is the 
single jump distance. For the diamond structure, A = 
0.25-\/3 (in units of the unit cell dimension) and so we 
obtain the correlation factor to be 



{xl) 



3 x (0.25) 2 x (N) 



(14) 



From the output of the KMC simulation, which has the 
displacements and the number of jumps of each atom, 
the mean squared displacements ((A 2 ), (A 2 ), (A 2 )) and 
the mean number of jumps ((A)) can be calculated by 
averaging the quantities over all atoms of the same type 
(Si or Ge) . The correlation factor can thus be calculated 
in a straightforward manner from the KMC simulation 
results. 



III. RESULTS AND DISCUSSION 

A. Biaxial strain dependence of the activation 
energy 

The change in vacancy formation energy due to biaxial 
strain is computed asE3 (— 2/3)V r /i, where V r is the relax- 
ation volume and fj, is the biaxial modulus of Si. From our 
LDA calculations, the relaxation volume accompanying 
the formation of a vacancy in Si is -20.73 A 3 . (The re- 
laxation volume in this case is approximately -\~Q6 times 
the atomic volume of Si.) The biaxial modulusEa of Si is 



7 



> 
>- 
c 

UJ 



u 

I— 

<v 
c 




• 


% Strain 


▲ 


+ 0.4% Biaxial 


T 


-0.4% Biaxial 



112 3 Far(12) 

Coordination Site from Ge 

FIG. 4: Interaction potential energy (in eV) between a sub- 
stitutional Ge atom and a vacancy as a function of vacancy 
position in (a) relaxed Si (circles connected with solid line) (b) 
0.4% tensile biaxially strained Si (upward triangles connected 
with dashed line) (c) 0.4% compressively biaxially strained Si 
(downward triangles connected with dotted line) from LDA 
calculations. The energies of all the systems when the va- 
cancy and the Ge are far apart have been normalized to zero. 
(Lines are drawn as a guide to the eye.) 



190.48 GPa. We therefore find the change in vacancy for- 
mation energy due to equi-biaxial strain to be 16 eV/unit 
strain. 

The interaction potential energies of a vacancy with a 
substitutional Ge atom for (a) relaxed, (b) 0.4% tensile 
equi-biaxially strained, and (c) 0.4% compressively equi- 
biaxially strained systems are as shown in Fig. |I| (The 
interaction potential energy for the relaxed system shown 
in Fig. |^ is exactly same as that shown in Fig. ||. It has 
been shown separately in Fig. || for the sake of clarity and 
has been shown in Fig. ^ for the sake of making compar- 
isons with the strained systems.) We make the follow- 
ing comments and observations with reference to Figs. ^ 
and |]: (i) Our first principles calculations indicate that 
a biaxial tension (compression) of 0.4% causes the third 
dimension to contract (expand) by 0.39% (0.46%) in rea- 
sonabljg-agreement with that predicted by linear elasticity 
theoryEa which gives a contraction (expansion) of 0.31% 
(0.31%). To maintain consistency, the interaction poten- 
tial energy calculations were made using the dimensions 
obtained from our LDA computations, (ii) From Fig. ^ 
we see that the migration barrier for the vacancy to ex- 
change positions with a Si atom far away from a Ge atom 
is 0.11 eV. This is as expected, owing to the similarity of 
the calculation technique, in reasonable agreement with 
Nelson et alE3 who report a value of 0.18 eV from their 
LDA calculations, (iii) The asymmetric location of the 
saddle points between the 1st and the 2nd coordination 



sites, and the 2nd and the 3rd coordination sites is due to 
the weaker nature of the Si-Ge bond compared to that of 
the Si-Si bond, (iv) From Fig. || we see that the binding 
energy of the vacancy to the Ge atom with the vacancy 
being at the n th coordination site from Ge is 0.24 eV, 
0.04 eV, and less than 0.002 eV respectively for n = 1, 2 
and 3. Therefore, the vacancy is practically bound to Ge 
only if it is at a nearest neighbor site to Ge. The strength 
and the range of interaction between a vacancy and Ge is 
quite weak unlike those between a vacancy and a dopant 
atom such as arsenicE3 or phosphorousCJ. This difference 
in the intensity and the extent of the interaction and the 
difference in the typical concentration of Ge in Sii-^Ge^ 
alloys compared to dopant concentrations suggests that 
the diffusion of Ge will not be dominated by . th e pair dif- 
fusion mechanism, which is the acceptedrnr ! dominant 
mechanism of dopant diffusion diffusing by the vacancy 
mechanism. Rather, the vacancy by randomly moving 
through the crystal randomly displaces Ge atoms when- 
ever it meets one, thereby causing diffusion. It does not 
form as strong a pair with the Ge atom as, for exam- 
ple, it does with a phosphorous or an arsenic atom, (v) 
From the interaction potentials (Fig. |J), we find that the 
barrier for the Si-V jump (far from a Ge atom) changes 
by 4eV/unit equi-biaxial strain and the barrier for the 
Ge-V jump (at very low Ge concentration) changes by 
2eV/unit equi-biaxial strain. 

The Ge-V interaction potential from Fig. ^ can be used 
to calculate the correlation factor for Ge diffusion as 
outlined in Sec 
(TST) 



[ID 



From the transition state theory 
, the jump rates i>i, uh, vp and v B mentioned in 



Sec. II D can be calculated as: 



vi = v exp[-(E xs ~ Ei)/k B T] 

vh = foexp[-(I?/ s - Ef)/k B T] 

v F = j/ exp[-(£'i2 - E 2 )/k B T] 

v B = ^oCxp[-(£i2 - Ei)/k B T] 



(15) 
(16) 
(17) 
(18) 



where vq denotes the lattice vibrational frequency (which, 
to a first order approximation we have assumed to be a 
constant); E\, E% and Ef respectively denote the energy 
of the system when the vacancy is at the first nearest 
neighbor site to the Ge atom, second nearest neighbor 
site to the Ge atom and far away from the Ge atom; E xs , 
Ef s and E\2 respectively denote saddle point energies 
for the vacancy and Ge to exchange positions, for the 
vacancy and Si to exchange positions far away from a Ge 
atom, and for the vacancy to move between the first and 
the second nearest neighboring positions of the Ge atom. 
Fig. [| plots the variation of the correlation factor for 
Ge diffusion as aJpnction of temperature. (As we have 
noted previouslyEJ the correlation factor approaches the 
theoretical limit of 0.5 at high temperatures.) In Fig. ^ 
we show an Arrhcnius plot of the same and we extract the 
activation energy for the strain free, 0.4% biaxial tensile 
and 0.4% biaxial compressive cases to be 0.168 eV, 0.171 
eV, and 0.161 eV respectively. The activation energy 
associated with the correlation factor for Ge diffusion in 
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FIG. 5: Theoretical calculation of the correlation factor for 
the diffusion of Ge in (a) relaxed Si (solid line) (b) 0.4% tensile 
biaxially strained Si (dashed line) (c) 0.4% compressively bi- 
axially strained Si (dotted line) as a function of temperature. 
The correlation factors are seen to approach a high tempera- 
ture limit of 0.5, the theoretical value for a tracer diffusion in 
a diamond structure. 
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FIG. 6: An Arrhenius type plot of the correlation factor for 
the diffusion of Ge in (a) relaxed Si (solid line, E c ° rr = 0.168 
eV) (b) 0.4% tensile biaxially strained Si (dashed line, E c ° rr = 
0.171 eV) (c) 0.4% compressively biaxially strained Si (dot- 
ted line, Ea° rr = 0.161 eV) to extract the activation energy 
corresponding to the correlation factor (E^ orr ). 



Si at very low Ge concentrations therefore changes by 
approximately 1 eV/unit strain. 

Combining the results of the vacancy formation en- 
ergy change due to biaxial strain with the migration bar- 
rier energy change and the correlation factor activation 
energy change, we estimate the following values for the 
effect of equi-biaxial strain on the diffusion-activation en- 
ergy: 20 eV/unit strain for Si self diffusion in Si, 17 - 20 



eV/unit strain for Ge self diffusion in Si. At this point, 
we would like to quote two experimentally determined 
values and one empirically fitted value for the change in 
activation energy for Ge diffusion due to biaxial strain. 
The experimentally determined, values are due to Cowern 
et alE3 and Zangenberg et al£3 who report a value of 18 
eV/unit strain and 160 ± 40 eV/unit strain respectively. 
Aubertine et aln, who use the strain dependence of the 
activation energy as a tunable parameter in their empir- 
ical model, report that they are best able to reproduce 
their experimental data if they set this parameter to be 
19eV/unit strain. 



B. Ge concentration dependence of the vacancy 
formation energy 

The typical concentration of Ge in SiGe films in device 
structures is 10% - 30%, which is several orders of magni- 
tude larger than the typical dopant concentration (10 16 
- 10 18 per cm 3 ). One therefore needs to consider the ef- 
fect of Ge concentration on the vacancy formation energy 
(and hence the vacancy conc entra tion) in the system in 
the manner explained in Sec. II C. From straightforward 
LDA based calculations, we obtain the following energet- 
ics of the SiGe- vacancy system: The energy of the system 
drops by 0.24, 0.45, 0.6, and 0.83 eV when the vacancy is 
surrounded respectively by 1, 2, 3 and 4 Qe. atoms. Simi- 
lar figures have been previously reported.OLJ The energy 
of the system drops by 0.04 eV for every second nearest 
neighbor Ge atom to the vacancy. The energy of the sys- 
tem with Ge in all the second nearest neighbor positions 
of the vacancy and Si everywhere else is higher than the 
energy of a system with a vacancy in un ary G e by 0.12 
eV. In terms of the notations used in Sec. IIC, F(0) = 
eV, F(l) 



0.24 eV, F{2) = 0.45 eV, F(3) = 0.6 eV, 
F(4) = 0.83 eV, S = 0.04 eV, and M = 0.12 eV. The 
attractive interaction between the Ge atoms and the va- 
cancy causes the equilibrium vacancy concentration to be 
larger in regions of high Ge concentration. This lowers 
the effective vacancy formation energy in Sii-^Ge^ com- 
pared to a uniform (regionally unbiased) random distri- 
bution of vacancies. 

Using the theory explained in Sec. [II C| , we have plotted 
in Fig. [7] the change in the vacancy formation energy from 
the vacancy formation energy in Si ((Ef(x,T)) — E v l f ) 
as a function of Ge concentration calculated at 1000K 
(solid line). Also plotted in Fig. [j] is the change in the 
change in the vacancy formation energy vs Ge concentra- 
tion from a rule-of-mixtures model for the composition 
dependence of the formation energy (dotted line). The 
rule-of-mixtures model is consistent with a spatially uni- 
form distribution of vacancies for each Sii-^Ge^ com- 
position. The difference between the two curves has a 
maximum at a particular concentration of Ge, which, of 
course, is temperature dependent. This is understood by 
the following reasoning: At very high Ge concentrations, 
a randomly chosen site would have a high probability 
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Ge Concentration 

FIG. 7: Solid line shows the change in the vacancy for- 
mation energy in Sii-^Ge^ from that in pure Si (dEf = 
(Ef(x,T)) — EyV) as a function of Ge concentration (x) cal- 
culated at 1000K by taking into account the attractive inter- 
action of the vacancy with the Ge atoms. Dotted line shows 
the same quantity obtianed by a simple rule of mixtures. 

of having many Ge neighbors. The further reduction in 
the formation energy because of the vacancies preferen- 
tially forming at high Ge concentration sites is therefore 
marginal. At very low Ge concentrations, the amount of 
reduction in the formation energy is low because of the 
small number of Ge atoms present. 

C. Vacancy migration energy barrier database 

We present, in this subsection, the database of energy 
barriers for vacancy migration in different environments 
calculated using the local density approximation. As in 
the case of calculating the effective vacancy formation en- 
ergy, we have treated atoms at different distances from 
the vacancy migration center differently depending on 
the extent of influence that the atom would exert on the 
vacancy migration energy barrier. Referring to Fig. ||, 
the identities of the atoms that are first nearest neigh- 
bors to the vacancy (denoted as SI, S2, S3 in the fig- 
ure), the identity of the migrating atom (denoted as DO), 
and the identities of the atoms surrounding the migrat- 
ing atom (denoted as Dl, D2, D3) are expected to have 
the greatest influence on the migration energy barrier. 
We have assumed that the concentration of vacancies is 
sufficiently low that none of the seven atoms (SI - S3, 
DO - D3) would be a vacancy. We then get a list of 40 
different configurations depending on which of (SI - S3, 
DO - D3) is Si or Ge. We account for the effect of the 
identities of the atoms beyond these seven nearest neigh- 
bors in the following mean field manner: We calculate 
the migration energy barriers for the 40 different config- 
urations for the following two cases: (i) All the atoms 
beyond the seven nearest neighbors are Si. (ii) All the 
atoms beyond the seven nearest neighbors are Ge. Then, 




Initial Final 



FIG. 8: The energy barrier for the vacancy (filled square) to to 
go from the initial configuration to the final is influenced the 
most by the identities of the atoms surrounding the vacancy 
(SI, S2, S3), the identity of the atom with which the vacancy 
is to exchange position (DO), and the identities of the atoms 
surrounding DO, namely, Dl, D2, and D3. A two dimensional 
representation of the diamond structure has been adopted 
for convenience with the different types of lines representing 
bonds on different planes. 



to obtain the energy barrier for any one of these 40 con- 
figurations in a Sii-aGe^ alloy (with a Ge concentration 
of x), we linearly interpolate the migration energy bar- 
rier of that particular configuration from cases (i) and 
(ii) mentioned above. 

This approach seems to be reasonably satisfactory for 
at least one of the configurations (all (SI - S3, DO - D3) 
are Ge atoms) that we have tested (Fig. ^ . The top con- 
figuration shown in Fig. ^ corresponds to the case where 
the seven nearest neighbors within the dotted circle (SI - 
S3, DO - D3) are all Ge atoms and all atoms beyond the 
nearest neighbor sites are Si (case (i) above) . The migra- 
tion energy barrier is 0.03 eV. The bottom configuration 
corresponds to the same nearest neighbor configuration 
but all atoms beyond the nearest neighbor sites being Ge 
(case (ii) above). The migration energy barrier is 0.13 
eV. The middle configuration is an explicit calculation of 
the energy barrier with the same nearest neighbor con- 
figuration (i.e., all (SI - S3, DO - D3) are Ge atoms) but 
all the atoms beyond the nearest neighbor sites are ei- 
ther Si or Ge with a probability of 0.5. This is consistent 
with what would occur in a Sii-^Gea; alloy with x = 0.5. 
Unlike in the top and the bottom configurations, the bar- 
rier for the forward migration (0.11 eV) is different from 
that for the reverse migration (0.07 eV) because unlike 
in the top and the bottom configurations, the identities 
of the second nearest neighbor sites are not identical for 
the forward and the reverse migrations. The mean of 
the forward and the reverse barriers (0.09 eV), however, 
is closer to the linear interpolation of the barriers from 
the top and the bottom configurations (0.08 eV) than it 
is to cither of the them. We should mention that the 
Sio.sGeo.s case is one which is farthest away from the 
reference cases (cases (i) and (ii)) and so is a stringent 
test, in a certain sense, of the approximation used. So 
a reasonable agreement in this case definitely indicates 
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FIG. 9: Although all the three configurations shown above 
have the same nearest neighbor atoms to the vacancy migra- 
tion center (atoms within the dotted circle), they have differ- 
ent vacancy migration energy barriers. The top configuration 
has Si atoms (open circles) everywhere outside the dotted cir- 
cle and has a vacancy migration energy barrier of 0.03 eV. 
The bottom configuration has Ge atoms (filled circles) every- 
where outside the dotted circle and has a vacancy migration 
energy barrier of 0.13 eV. The middle configuration has 50% 
of the atoms outside the dotted circle as Si. The barrier for 
the forward migration is 0.11 eV and that for the reverse mi- 
gration is 0.07 eV. The mean (0.09 eV) is closer to the linear 
interpolation of the barriers from the top and the bottom con- 
figurations (0.08 eV) than it is to either of them. (Si: open 
circles; Ge: filled circles; vacancy: filled square with a label 
"V".) 



that the approximation used is reasonable. 

We have calculated the saddle point energies for each 
of these 80 different configurations very accurately using 
the nudged elastic band (NEB) method.EJ We also note 
that we have eliminated strain effects by suitably adjust- 
ing the lattice parameter for the number of Si and Ge 
atoms for each of the 80 configurations. Table | summa- 
rizes the vacancy migration energy barriers. The barriers 
are negligibly small for some of the configurations where 
there is a significant asymmetry between the initial and 
the final environments of the vacancy (in terms of the 
number of Ge atoms) especially for the case of 0% Ge. 
By plotting the atomic positions, we have found out that 
the reason for this behavior is the following: When there 
is a significant assymetry, the initial (or the final) struc- 
ture "collapses" to the other or to some configration in- 
termediate between the two and there is no barrier to get 
from the one to the other. We feel that it happens more 
in the case of 0% Ge because the lower lattice constant of 
Si compared to Ge facilitates this collapsing more easily 
than in the case of 100% Ge. 



D. Kinetic Monte Carlo simulation 

In this section we present the results of the kinetic 
Monte Carlo simulations done using the vacancy migra- 
tion energy barrier database presented in Sec. Ill C . 



1. Diffusivity, AMF energy 

From the KMC simulations, we were able to compute 
the diffusivity of Si and Ge in Sii-zCe^ as a function 
of Ge concentration (x) and the temperature. (The dif- 
fusivity is given by D = (X 2 )/6t, see Sec. [II E| ). We 
note that these simulations have a constant vacancy con- 
centration (of 10 _6 /atom); in other words, the change 
in the vacancy concentration due to the c hange in the 
Ge concentration as explained in Sec. Ill B has not been 
factored in. The lattice vibrational frequency vq was esti- 
mated from first principles based on a harmonic approxi- 
mation to be 7.325 x 10 11 sec -1 . Figures [To|and [PI show 
the plots of the Ge and the Si diffusivities respectively. 
As expected, the diffusivity increases with temperature. 
From an Arrhenius type plot of the diffusivities, we have 
extracted the activation minus formation (AMF) energy. 
These have been plotted in Fig. [l2| We make the fol- 
lowing observations with reference to Figs. [Io| - [TJ. (i) 
The AMF energy for the diffusion of Si in Si (0.11 eV) 
matches closely with the migration energy for a vacancy 
in pure Si. (Compare with the entry in Table [I] corre- 
sponding to all (SI - S3, DO - D3) being Si under 0% 
Ge.) A similar close match is also obtained for the AMF 
energy for the diffusion of Ge in Ge (0.13 eV). (Compare 
with the entry in table 1 corresponding to all (SI - S3, 
DO - D3) being Ge under 100% Ge.) Along with provid- 
ing a verification of our computer simulation programs, 
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TABLE I: Vacancy migration energy barrier database. SI, 
S2, and S3 are the identities of the atoms surrounding the 
vacancy, DO is the identity of the atom with which the vacancy 
is to exchange positions, and Dl, D2, and D3 are the identities 
of the atoms surrounding DO. Under 0% Ge are listed the 
energy barriers corresponding to the case where all atoms 
other than SI - S3, DO - D3 are Si and under 100% Ge are 
listed the barriers when those atoms are all Ge instead. 



SI 


S2 


S3 


DO 


Dl 


D2 


D3 


0% Ge (eV) 


100% Ge (eV) 


Si 


Si 


Si 


Si 


Si 


Si 


Si 


0.11 


0.26 


Si 


Si 


Si 


Ge 


Si 


Si 


Si 


0.08 


0.22 


Si 


Si 


Si 


Si 


Si 


Si 


Ge 


0.03 


0.87 


Si 


Si 


Si 


Ge 


Si 


Si 


Ge 


0.00 


0.10 


Si 


Si 


Si 


Si 


Ge 


Ge 


Si 


0.00 


0.05 


Si 


Si 


Si 


Ge 


Ge 


Ge 


Si 


0.00 


0.03 


Si 


Si 


Si 


Si 


Ge 


Ge 


Ge 


0.00 


0.00 


Si 


Si 


Si 


Ge 


Ge 


Ge 


Ge 


0.00 


0.00 


Si 


Si 


Ge 


Si 


Si 


Si 


Si 


0.23 


1.03 


Si 


Si 


Go 


Ge 


Si 


Si 


Si 


0.18 


0.81 


Si 


Si 


Go 


Si 


Si 


Si 


Ge 


0.09 


0.86 


Si 


Si 


Ge 


Si 


Ge 


Si 


Si 


0.10 


0.89 


Si 


Si 


Ge 


Ge 


Si 


Si 


Ge 


0.06 


0.18 


Si 


Si 


Ge 


Ge 


Ge 


Si 


Si 


0.06 


0.67 


Si 


Si 


Ge 


Si 


Ge 


Ge 


Si 


0.01 


0.74 


Si 


Si 


Ge 


Si 


Si 


Ge 


Ge 


0.00 


0.70 


Si 


Si 


Ge 


Ge 


Ge 


Ge 


Si 


0.00 


0.09 


Si 


Si 


Ge 


Ge 


Si 


Ge 


Ge 


0.00 


0.07 


Si 


Si 


Ge 


Si 


Go 


Ge 


Ge 


0.00 


0.64 


Si 


Si 


Go 


Ge 


Go 


Ge 


Ge 


0.00 


0.01 


Si 


Go 


Ge 


Si 


Si 


Si 


Si 


0.00 


0.51 


Si 


Ge 


Ge 


Ge 


Si 


Si 


Si 


0.00 


0.44 


Ge 


Ge 


Si 


Si 


Si 


Si 


Ge 


0.21 


0.34 


Ge 


Ge 


Si 


Si 


Ge 


Si 


Si 


0.00 


0.32 


Ge 


Ge 


Si 


Ge 


Si 


Si 


Ge 


0.16 


0.29 


Ge 


Ge 


Si 


Ge 


Ge 


Si 


Si 


0.00 


0.27 


Go 


Go 


Si 


Si 


Ge 


Ge 


Si 


0.07 


0.18 


Go 


Go 


Si 


Si 


Si 


Ge 


Ge 


0.08 


0.20 


Ge 


Go 


Si 


Ge 


Go 


Ge 


Si 


0.04 


0.15 


Ge 


Go 


Si 


Ge 


Si 


Ge 


Ge 


0.05 


0.16 


Si 


Go 


Ge 


Si 


Go 


Ge 


Ge 


0.00 


0.08 


Si 


Ge 


Ge 


Ge 


Ge 


Ge 


Ge 


0.00 


0.05 


Ge 


Ge 


Ge 


Si 


Si 


Si 


Si 


0.00 


0.69 


Ge 


Ge 


Ge 


Ge 


Si 


Si 


Si 


0.00 


0.00 


Ge 


Ge 


Ge 


Si 


Si 


Si 


Ge 


0.00 


0.48 


Go 


Go 


Ge 


Ge 


Si 


Si 


Ge 


0.00 


0.41 


Go 


Go 


Ge 


Si 


Si 


Ge 


Ge 


0.00 


0.30 


Ge 


Ge 


Ge 


Ge 


Si 


Ge 


Ge 


0.00 


0.25 


Ge 


Go 


Ge 


Si 


Go 


Ge 


Ge 


0.05 


0.16 


Ge 


Go 


Go 


Ge 


Go 


Ge 


Ge 


0.03 


0.13 



it also corroborates our concept of the AMF energy as 
explained in Sec. |. (ii) We have not been able to obtain 
satisfactory numerical agreement of the AMF energy for 
the diffusion of Si in Ge (0.34 eV) or for that of diffu- 
sion of Ge in Si (0.05 eV) based on the models presented 
in Refs. plp^p3| . We however have the following plau- 
sible qualitative explanation: The attractive interaction 
between vacancies and Ge atoms (see Fig. ||) causes a va- 
cancy to be more available near the vicinity of a Ge atom 
to facilitate its diffusion. This probably results in lower- 
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FIG. 10: Diffusivity of Ge in Sii-zCe^, calculated from the 
results of the KMC simulation as a function of Ge concentra- 
tion (x) at five different temperatures: 300K - square; 600K 
- circle; 900K - upward triangle; 1200K - downward triangle; 
1500K - pentagram. (Lines are drawn as a guide to the eye.) 



ing the AMF energy for the Ge diffusion in Si compared 
to the migration barrier for a Ge- vacancy exchange pro- 
cess in Si (0.08 eV (see the second entry under 0% Ge in 
Table |)). Conversely, the repulsive interaction between 
a Si and a vacancy (see Fig. |l3|) makes a vacancy less 
available near the vicinity of a Si atom. This probably 
results in increasing the AMF energy for the Si diffusion 
in Ge compared to the migration barrier for a Si- vacancy 
exchange process in Ge (0.16 eV (see the penultimate en- 
try under 100% Ge in Table |)). (iii) While we do not 
have a microscopic explanation for the abrupt drop in the 
diffusivity of both Si and Ge near low Ge concentrations, 
we do find them to be consistent with the rise in AMF 
energy of both Si and Ge near low Ge concentrations, 
(iv) The reason for the non smooth behavior of the AMF 
energies near 50% Ge concentration is probably because 
of those concentrations being farthest away from the ref- 
erence configurations (0% and 100% Ge) that were used 
to build the migration energy barrier database. 

In Fig. |lj, we plot the change in the activation energy 
for Ge diffusion in Sii_ x Ge x compared to that in Si as 
a function of Ge concentration. The activation energy is 
calculated as a sum of the vacancy formation energy from 
Fig. ^ and the Ge AMF energy from Fig. Also plot- 
ted on the same axes are the experimentally observed 
changes in the activation energy for Ge diffusion from 
Refs. [Tl],[l2|,|l3]. The purpose of this plot is not actually 
to compare our results with the experiments, which, as 
we have already mentioned in Sec. |j is premature at this 
stage. But this plot clearly brings out the need to con- 
sider the other mechanisms before a fair comparison with 
experiments is possible. 
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FIG. 11: Diffusivity of Si in Sii-^Ge^, calculated from the 
results of the KMC simulation as a function of Ge concentra- 
tion (x) at five different temperatures: 300K - square; 600K 
- circle; 900K - upward triangle; 1200K - downward triangle; 
1500K - pentagram. (Lines are drawn as a guide to the eye.) 
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FIG. 13: Interaction potential energy (in eV) between a sub- 
stitutional Si atom and a vacancy as a function of vacancy 
position in relaxed Ge from LDA calculations. (Lines are 
drawn as a guide to the eye.) 
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FIG. 12: Variation of the activation-minus-formation (AMF) 
energy (eV) for the diffusion of Ge (upward triangle) and Si 
(circle) in Sii-^Ge^ as a function of Ge concentration (x) 
obtained from the results of the KMC simulation. (Lines are 
drawn as a guide to the eye.) 
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FIG. 14: Solid line shows the change in the activation energy 
for the diffusion of Ge in Sii-zCe^ by a vacancy mechanism 
from that in pure Si as a function of Ge concentration (x) cal- 
culated at 1000K as the sum of the change in the vacancy for- 
mation energy and the Ge activation-minus-formation (AMF) 
energy. Also shown are the expcppiental results for the same 
quantity from Zangenberg et al!El (upward itaiangle) , Strohm 
et ali£3 (circle) , and McVay and DuCharmetj (downward tri- 
angle) . 



The correlation factor for Ge a nd Si cal culated by the 
procedures outlined in Sees. II D| and [IE are plotted re- 
spectively in Figs. [[5] and |16| as a function of the Ge 
concentration for five different temperatures. We make 
the following observations with reference to these plots, 
(i) The correlation factors that we have calculated for 
the unary substances (i.e., Si correlation factor in 0% 



Ge concentration and Ge correlation factor in 100% Ge 
concentration) equal 0.5. This is the theoretical value 
for the correlatioipJactor for a tracer diffusion in the di- 
amond structure. E3 (This provides an additional verifi- 
cation of our calculations.) That the correlation factor 
for a tracer (diffusing by the vacancy mechanism) should 
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be less than unity is clear from the observation that the 
tracer atom has a higher probability of jumping back to 
the vacancy site thereby nullifying a forward jump. The 
mean squared displacement (and hence the diffusivity) of 
this correlated motion would therefore be less than that 
of a random jump, giving a correlation factor less than 
unity. The specific value of 0.5 is a result of the tetra- 
hedral geometry of the silicon crystal structure, (ii) At 
higher temperatures, the Boltzmann factor evens out the 
different energy barriers, making the system resemble a 
unary substance. One would therefore expect the corre- 
lation factor to approach the value for a unary substance 
in the diamond structure, namely 0.5. We observe this 
in our plots, (iii) At low Ge concentrations, the correla- 
tion factor for Ge drops below 0.5. The reason for this 
is understood by the following argument: The attractive 
interaction between Ge and a vacancy and the lower en- 
ergy barrier for a vacancy to exchange positions with Ge 
than to jump to the second nearest neighbor site of Ge 
from the first (which results in breaking the Ge-vacancy 
bond, see Fig. ||), tend to cause the Ge and the vacancy to 
jump back and forth several times before breaking away 
from each other. But, in the diamond structure, because 
there is no atomic location that is a simultaneous neigh- 
bor to both the vacancy and the Ge atom (when they 
are first nearest neighbors to each other), breakage of the 
Ge-vacancy pair is essential for the Ge atom to be effec- 
tively displaced from its current location. This back and 
forth motion does not contribute to the mean squared 
displacement of the Ge atoms and consequently the Ge 
correlation factor drops, (iv) At low Ge concentrations, 
the correlation factor for Si drops below 0.5. We offer 
the following explanation for this behavior: The attrac- 
tive interaction between the Ge and the vacancy causes 
the vacancies to be predominantly found near Ge atoms. 
So, only the Si atoms found near those Ge atoms are af- 
fected by the vacancy motion. These Si atoms, owing to 
the lower energy barrier for the vacancy to jump to the 
first nearest neighbor site of the Ge atom from the second 
than to jump to the third from the second (see Fig. ||), 
just keep jumping back and forth between the first and 
the second nearest neighbor sites of the Ge (depending on 
whether the vacancy is correspondingly at the second or 
the first nearest neighbor sites) . This back and forth mo- 
tion does not effectively displace the Si atoms and so does 
not contribute to the mean squared displacement. This 
causes the Si correlation factor to drop, (v) The corre- 
lation factors of both Si and Ge increase with increasing 
Ge concentration. We explain this in the following man- 
ner: As the concentration of Ge increases, the vacancy 
is attracted by the other Ge atoms too and therefore it 
is less likely to be bound to a single Ge atom. This re- 
duces the redundant back and forth motion of the Ge 
atoms, thus increasing the mean squared displacement 
and consequently the Ge correlation factor. (This effect 
is similar, iiL-spme ways, to the percolation mechanism 
for diffusion .E!l) The vacancy, in the process of moving 
from one Ge atom to the other, ends up displacing Si 
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FIG. 15: Correlation factor for the diffusion of Ge in Sii-^Ge^ 
calculated from the results of the KMC simulation as a func- 
tion of Ge concentration (a;) at five different temperatures: 
300K - square; 600K - circle; 900K - upward triangle; 1200K 
- downward triangle; 1500K - pentagram. (Lines are drawn 
as a guide to the eye.) 



atoms thus increasing their mean squared displacement 
and consequently the Si correlation factor, (vi) The cor- 
relation factor for Si in Sii-^Ge^ alloys with high Ge 
concentration is greater than 0.5 and approaches unity. 
This interesting behavior is explained by the following 
reasoning: At very high Ge concentration (i.e., very low 
Si concentration), the faster jumping rate of the vacancy 
with the Ge atoms compared to that with the Si atoms 
(because of the lower barrier height (compare last two en- 
tries under 100% Ge in Table |)) causes the vacancy to 
perform a lot of jumps with Ge atoms between successive 
jumps with a Si atom. This results in the vacancy ap- 
proaching Si via an essentially random path, making the 
Si jumps closer to a random walk process. This causes 
the correlation factor to approach unity. In Fig. [l3|, we 
show the interaction energy between a Si and a vacancy 
in a Ge environment. In Fig. [l?] we show the variation 
of the Si correlatio n fac tor with temperature calculated 
as outlined in Sec. 



II D 



We do see that the correlation 



factor tends to unity in the lower temperature limit. 



IV. SUMMARY 

Our purpose of the present work was to understand, 
from first principles, the effect of biaxial strain and com- 
position on the self-diffusivity of Si and Ge in Sii-^Ge^ 
alloys. In order to attack the problem, we broke it down 
into one of studying the effect of these factors on the 
main components that define the diffusivity: the va- 
cancy formation energy, and the activation minus forma- 
tion (AMF) energy. (The necessity and the definition of 
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FIG. 16: Correlation factor for the diffusion of Si in Sii-xGe^ 
calculated from the results of the KMC simulation as a func- 
tion of Ge concentration (x) at five different temperatures: 
300K - square; 600K - circle; 900K - upward triangle; 1200K 
- downward triangle; 1500K - pentagram. (Lines are drawn 
as a guide to the eye.) 
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FIG. 17: Theoretical calculation of the correlation factor for 
the diffusion of Si in relaxed Ge. 



AMF energy were presented.) We attacked the problem 
by the following three main steps: (i) We performed den- 
sity functional theory (DFT) calculations within the local 
density approximation (LDA) to obtain the required en- 
ergetics of the various configurations, (ii) We worked out 



the details necessary to calculate the correlation factor 
and the change in the vacancy formation energy with 
composition, (iii) We performed kinetic Monte Carlo 
(KMC) simulations using our total energy calculations. 
By this approach, we were able to estimate the following 
values for the effect of biaxial strain on the activation 
energy (the sum of the vacancy formation energy and 
AMF energy): 20 eV/unit strain for Si self diffusion in 
Si and 17 - 20 eV/unit strain for Ge self-diffusion in Si. 
We calculated the change in the vacancy formation en- 
ergy in Sii-zGe^ as a function of composition. From 
the KMC simulations, we were able to extract the vari- 
ation of the AMF energy for Si and Ge self-diffusion in 
Sii-zGea; as a function of composition. We combined 
the Ge AMF energy with the vacancy formation energy 
to find the variation of the activation energy for Ge dif- 
fusion in Sii-zGejc as a function of composition. Lastly, 
we presented the variation of the correlation factor for 
Si and Ge diffusion in Sii-^Ge^ as a function of com- 
position and temperature and made several interesting 
observations that are quite general for a vacancy medi- 
ated diffusion in a random binary alloy arranged in a 
diamond structure. 

There are many outstanding issues of the complete 
model that need to be resolved even for the vacancy 
mechanism alone. We conclude this article by recogniz- 
ing the following limitations of the present work: (i) As 
we mentioned in the introduction, the inability of the 
LDA to reproduce experimentally observed values of the 
activation energy in Si precludes our results from being 
directly compared with experiments, (ii) We have not 
addressed the effect of strain and composition on the 
pre-exponential factor and have not considered entropic 
effects. 
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